function out = forecastRegressionYieldSpreadPCA_forBoot(yields,shortRateForXhr,shortRateOneYear,maturities,REC,h)
% Data inputs:
% yields: T * numYields where yields are annually and continuously compounded
% shortRateForXhr: Short rate for excess returns
% maturities: 1*numYields with maturities in months
% RECdummy: 0 1 regression dummy, where 0 is boom and 1 is recression
% h: sampletime_per_maturityTime

% Preable
T         = size(yields,1);
numYields = length(maturities);

% The PCA's
[Wpca,Xpca] = pca(yields,'NumComponents',3,'Centered',false);
numPCAs     = size(Xpca,2);

% Term spread: y[t](n)-y[t](1-year) for all n.
termSpreads      = yields(:,1:end)-repmat(shortRateOneYear,1,numYields);

% Ex post holding period returns: hr[t]_n = n*y[t-h](n) - (n-h)*y[t](n-h)
annualMaturities = maturities/12;
holdPerReturn   =[nan(T-h,1), (repmat(annualMaturities(1,2:end),T-h,1).*yields(1:end-h,2:end)...
    -repmat(annualMaturities(1,1:end-1),T-h,1).*yields(1+h:end,1:end-1))];

% Ex post excess holding period returns: xhr[t]_n = (hr[t]_t - y[t-h](h))*1/h           
excessHoldPerReturn = (holdPerReturn/h - repmat(shortRateForXhr(1:end-h,1)/12,1,numYields));               
               
%% Runing the forecast regressions for all maturities and for the three predictors
coef_alfa      = NaN(1,numYields-1);
coef_PCA       = NaN(numPCAs,numYields-1);
coef_dalfa     = NaN(1,numYields-1);
coef_dbeta     = NaN(1,numYields-1);

se_alfa        = NaN(1,numYields-1);
se_PCA         = NaN(numPCAs,numYields-1);
se_dalfa       = NaN(1,numYields-1);
se_dbeta       = NaN(1,numYields-1);

tstat_alfa     = NaN(1,numYields-1);
tstat_PCA      = NaN(numPCAs,numYields-1);
tstat_dalfa    = NaN(1,numYields-1);
tstat_dbeta    = NaN(1,numYields-1);

R2regime         = NaN(1,numYields-1);
startMat         = find(maturities >= 2*12);

%for n=1:numYields-1
for n=startMat:numYields
    % excessHoldPerReturn regressed on term spread
     Y = excessHoldPerReturn(:,n);
          
    % Regime dependent: excessHoldPerReturn regressed on constant, PCA's, recessionIntercept and recssionSlope
    X1=[ones(size(Y,1),1)  Xpca(1:end-h,:) REC(1:end-h,1) termSpreads(1:end-h,n).*REC(1:end-h,1) ];
    resRegime            = nwest(Y,X1,2*h);
    coef_alfa(1,n)       = resRegime.beta(1);
    coef_PCA(:,n)        = resRegime.beta(2:4);
    coef_dalfa(1,n)      = resRegime.beta(5);
    coef_dbeta(1,n)      = resRegime.beta(6);
    
    tstat_alfa(1,n)      = resRegime.tstat(1);
    tstat_PCA(:,n)       = resRegime.tstat(2:4);
    tstat_dalfa(1,n)     = resRegime.tstat(5);
    tstat_dbeta(1,n)     = resRegime.tstat(6);
    
    se_alfa(1,n)         = coef_alfa(1,n)/tstat_alfa(1,n);
    se_PCA(:,n)          = coef_PCA(:,n)./tstat_PCA(:,n);
    se_dalfa(1,n)        = coef_dalfa(1,n)/tstat_dalfa(1,n);
    se_dbeta(1,n)        = coef_dbeta(1,n)/tstat_dbeta(1,n);
    
    R2regime(1,n)         = resRegime.rsqr;
end

%% Output
out.coef_alfa       = coef_alfa;
out.coef_PCA        = coef_PCA;
out.coef_dalfa      = coef_dalfa;
out.coef_dbeta      = coef_dbeta;
out.tstat_alfa      = tstat_alfa;
out.tstat_PCA       = tstat_PCA;
out.tstat_dalfa     = tstat_dalfa;
out.tstat_dbeta     = tstat_dbeta;
out.se_alfa         = se_alfa;
out.se_PCA          = se_PCA;
out.se_dalfa        = se_dalfa;
out.se_dbeta        = se_dbeta;
out.R2regime         = R2regime;
out.maturities       = maturities;


end

